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Abstract 

Paper in honour of Freeman Dyson on the occasion of his 80th birthday. Normal N-body systems 
relax to equilibrium distributions in which classical kinetic energy components are l/2fcT, but, when 
inter-particle forces are an inverse cubic repulsion together with a linear (simple harmonic) attraction, 
the system pulsates for ever. In spite of this pulsation in scale, r(t), other degrees of freedom relax 
to an ever-changing Maxwellian distribution. With a new time, r , defined so that r 2 d/dt = d/dr 
it is shown that the remaining degrees of freedom evolve with an unchanging reduced Hamiltonian. 
The distribution predicted by equilibrium statistical mechanics applied to the reduced Hamiltonian 
is an ever-pulsating Maxwellian in which the temperature pulsates like r~ 2 . 

Numerical simulation with 1000 particles demonstrate a rapid relaxation to this pulsating equi- 
librium. 

1 Introduction 

In classical mechanics when N bodies interact with forces derived from a potential 

V = J>„, (1) 

n 

where V n is of inverse nth power in |xi — Xj|, the Virial theorem reads 

\ ^ = 2T + 5^nV„ = 2B + 5^(n-2)V n . (2) 

n n 

Here 

/ = S ' mi|xi — x| 2 = Mr 2 = 1/2 S ' S ' M^ 1 minij jx, — Xj | 2 , (3) 

i i j 

and indices i,j run over the particles. We notice that the term with n = 2 disappears from the second 
sum in equation @. Furthermore if the simple harmonic term is V-2 = \w 2 J^i< S M~ 1 mimj\x.i — 
xj| 2 then it can be re-expressed as u) 2 I/2. We now specialise to the problem in which only V2 and 
V-2 are present, so that any two particles of separation Tij repel each other as r~ 3 and attract like 
Tij. We consider this special problem because equation J5J now reads 
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so / vibrates harmonically about the value E/uj 2 . If this excitation is present initially it will continue 
vibrating at the same amplitude for ever, despite the complication of the r~ 3 repulsions of the 
particles. Multiplying equation (4) by dl /dt and integrating, 

h^Lf = 2EI-u 2 I 2 -M 2 C 2 , (5) 

where the last term is the integration constant. Using / = Mr 2 this becomes 

l/2(r 2 + C 2 r~ 2 +u 2 r 2 )M = E, (6) 

which may be compared with the energy of a particle of mass M orbiting with specific angular 
momentum C under a central simple harmonic force. Evidently 

d 2 r/dt 2 = £ 2 r" 3 - <Jr. (7) 

To save writing unimportant details hereafter we take all the masses equal, so that = m = M/N . 
The equation of motion for Xi now reads 

md 2 Xi/dt 2 = — mtj 2 (xi — x) — <9V2/<9xi, (8) 

where the simple harmonic forces were combined using equation (3). We now employ rescaled variables 
defined by 

Xi = (xi - x)/r. (9) 

Then 

d 2 x ; d 2 r dr dXi d 2 Xi 2 -3 2 n v . -3 2 d ( 2 dX ; \ 

Introducing a new 'time' r by d/dr = r 2 d/dt, we notice that Ct is the azimuth of the particle in the 
imaginary orbit introduced under equation (6). The equations of motion become 

md 2 Xi/dr 2 = -m£ 2 Xi - dV 2 /dX h (11) 

where V2 = r 2 Vi- Since V% is homogeneous of degree —2 in the Xi, one merely replaces the Xi by Xi 
to make V2 from V%. The result does not depend on r explicitly. Equation (11) is thus an autonomous 
equation for the evolution of the reduced variables Xi, but as functions of r rather than t. 
From their definition (9) the Xi are constrained so that both 

V X; = and ^2 X ' ? = N - ( 12 ) 

Now the E in equation (1) is the energy relative to the centre of mass since the I is measured in that 
frame, see (3). Our energy equation is therefore 

E = i^m[d(xi-x)/d£] 2 + K-2 + V2 = ^M(r 2 + wV) + r~ 2 [i ^ m(dX ; /dr) 2 +Va]. (13) 

Eliminating r via (6) and multiplying by r 2 we obtain 

^MC 2 = m (dXi/dr) 2 + Va. (14) 

So the "energy" of the reduced variables in r-time is MC 2 /2. Had we directly integrated equation 
(11) to find this energy, it would not have been obvious that the integration constant was zero. 

We are now in a position to state our problem in statistical mechanics. Given that the Xi must 
satisfy the constraints (12), what is their statistical equilibrium and how does any such equilibrium 
translate back into the eternally pulsating variables Xi and Xi? We find the equilibrium in section 2. 
In section 3 we demonstrate by numerical experiment with 1000 particles that a system started well 
away from that pulsating equilibrium relaxes to the predicted ever-pulsating equilibrium. Finally in 
section 4 we remark on the solutions of the corresponding problem in quantum mechanics. 

The problem is exceptional in that the normal dissipation of the basic breathing mode via Violent 
Relaxation is exactly absent. Nevertheless the other modes of the system do dissipate so Violent 
Relaxation is not totally absent. Although the long range simple harmonic force seems very artificial 
the net effect is exactly that found by Newton for the gravitational force within a homogeneous 
body. Such bodies have inspired many delightful studies by the great mathematicians including one 
by Freeman Dyson J2|. The special case of our problem in which V 2 is zero is the N-body problem 
exactly solved by Newton in the Principia 0. The pulsating equilibrium idea arose in our earlier 
generalisations of his work to other extraordinary N-body problems that are exactly soluble in both 
classical and quantum mechanics The special case when the harmonic force is absent and the 

particles are on a line is the exactly soluble Calogero model 
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2 Statistical Mechanics of the Reduced System 



The constraint Xi 2 = N tells us that the 3iV coordinates lie on a hypersphere of radius y/~N in 3N 
dimensions. Each of the three constraints Xi = gives a hyperplane through the hypersphere's 
centre so the first of them reduces the 3iV-sphere to a (3iV — l)-sphere; imposing the others as well 
reduces that to a (SN — 3)-sphere. If V2 were zero the representative point would move freely on such 
a hypersphere and with V2 present the motion is still a Hamiltonian one in the angles. The reduced 
system (11) still has a conserved total angular momentum. Writing X- for dXi/dr the conservation 
law is 

Xi x mXj = ^J(xi/r) x mr 2 — (xi/r) = xi x mi] = J. (15) 

Incorporating this constraint too, the reduced system has 3N — 7 degrees of freedom. For details of the 
angles on the hypersphere and their corresponding momenta see the appendix to [5]. Using Lagrange 
multipliers m/3*/2,7, S to impose the constraints of constant 'energy' MC 2 /2, J, and ^ Xi 2 ,the 
equilibrium distribution function factorises into a distribution of momenta mX' given by, 

/ocexp[-/3*mX' 2 /2-7.(Xx mX.% (16) 

and one dependent on the spatial coordinates. 

Writing 7 = — /3* S~2 and omitting a further function of position 

/aexp[-/3*m(X'-nxX) 2 /2]. (17) 

Returning to our original variables x = rX and writing x = v, we have X' = r[v — (r/r)x] where 
r(t) is the ever pulsating scale. In terms of our old variables / takes the form 

/ oc exp[-/3*mr 2 (v- u) 2 /2], (18) 

where u = (r/r)x — f2r~ 2 x x. We see that v is distributed Maxwellianly relative to the mean 
u(x,t).This mean moves with a time-dependent 'Hubble' flow superposed on a time-dependent rota- 
tion flr~ 2 . Furthermore the temperature of the distribution is time-dependent with f3*r 2 taking the 
place of the normal j3 so that the temperature is proportional to r~ 2 . We intentionally omitted the 
spatial distribution from the above as it contains all the complications of the problem. With both 
attractive and repulsive forces present we expect phase transitions of solids to liquids and gases even 
in classical physics without quantum phenomena. At low temperatures we expect a solid lattice but 
it can not be perfectly regular as the spacing must increase outward as the pressure decreases. None 
of this prevents the whole body undergoing large amplitude rescaling pulsations with the associated 
time-dependent rotation predicted above. At high enough temperatures V2 is always small compared 
with the kinetic energy and it can be neglected except as the means by which the system relaxes to 
its pulsating equilibrium whose density is then given by 

P = Po exp-[6Z 2 +q(X 2 + Y 2 )], (19) 

where X = (X,Y,Z). The coeficient q is best expressed in terms of a reduced omega u> 2 — m/3*Q, 2 
and takes the form 4q = 3 — ui 2 + -y/(3 — uj 2 ) 2 + 2a; 2 . The Lagrange multiplier 5 is S = q + ui 2 /2. 
Notice that 5 — » q — » 3/2 as Q. — > 0. From equation (19) we may calculate the moment of inertia 
and thence the total angular momentum is J = SIM/q. The other Lagrange multiplier is determined 
from m(3* C 2 = 3 + ui 2 /q. 



3 Simulating the Approach to the Pulsating Equilibrium 

Numerical simulations were carried out on a system of 1000 particles of equal mass with pair inter- 
actions depending on their separation nj, 

V{r ij ) = ^[rl+r^\, (20) 

using the method of molecular dynamics [7j. The unit of mass was defined so that the total mass 
M = 1000m was equal to one, and the equations of motion were integrated using the Verlet velocity 
algorithm with a time step of 0.001. With this choice of parameters the period of the oscillation of 
r 2 is 7r time units. The starting configuration was constructed from a cubic array of particles with an 
initial interparticle spacing of 0.1 and the origin of the coordinates was defined as the cube centre. 
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Velocities were chosen randomly from a Gaussian distribution, and adjusted so that the velocity of 
the centre of mass was equal to zero and the total energy E was equal to some specified value. The 
coordinate system was rotated so that the angular momentum was along the z direction. In order 
to study the approach to equilibrium, the initial velocity and spatial distributions were perturbed 
by scaling the z velocities and z positions by a factor of two. Such a scaling leaves the angular 
momentum along the z-axis and unchanged, 

The perpetual pulsating is illustrated in figure The top curve in this figure shows r 2 as a 
function of time measured in periods of n time units for periods 50-60 since the beginning of the 
simulation. The harmonic pulsation of r 2 is clear and both amplitude and phase are the same as 
at the beginning of the simulation. The lower two curves show the total potential energy and the 
contribution from V2, the inverse square term. It can be seen that the latter term is only important 
during the phase of the pulsation when the system is compressed so that r 2 is small. 

The relaxation toward equilibrium of both the shape of the cluster and the distribution of the 
peculiar velocities (v — u) = v p towards equilibrium are shown in figures [5] and [3] Figure ||5J shows 
the shape of the cluster as measured by the ratio z 2 / ^ x 2 ) 1 ^ 2 , which tends to one at equilibrium, 
since the only angular momentum in the simulation is the small one statistically generated in the 
initial conditions. Although the system is initially far from equilibrium, the shape relaxes over 
about 5 periods to an approximately spherical distribution with equal second moments in x and z, 
although both quantities are pulsating. Fluctuations in the ratio remain for many periods. Figure 
© shows the changes in r 2 ^ Vp Z and r 2 ^ Vp X as functions of time. According to equation (18) 
these quantities should be constant and equal at equilibrium. Relaxation of the difference occurs over 
about 5 pulsation periods. 

The Maxwellian distribution of peculiar velocities (v p ) predicted in equation (18) is illustrated 
in figure The top part of the figure shows the velocity distributions at eight different phases of 
the pulsation. At the phase when the cluster is compact the velocity distribution is broader than 
when the cluster is extended, but in every case the distribution is Maxwellian. The middle part of 
the figure verifies that the the width of the distribution is inversely proportional to r as predicted 
by equation (18), as it shows that, when the velocities are rescaled by r, the distributions coincide. 
To decrease numerical fluctuations all these graphs average together the distributions of the x,y and 
z components of velocity relative to the predicted mean and also average the distributions taken at 
the same phases of fifteen pulsations at the end of the run. The logarithmic plot of the distribution 
function against (v px ) 2 is shown to be linear in the final graph, confirming the Maxwellian form of 
the distribution. 

4 Conclusions 

The predicted pulsating equilibrium is approached and forms a limit cycle. Limit cycles are well- 
known in dynamics but it is more unusual to come across a pulsating Maxwellian distribution that 
continues with no change of entropy. However a related case is found in the Planck distribution of 
photons which remains rescaled but unchanged in shape or entropy as the universe expands. That is 
not true of the distribution function of massive particles which need collisions to maintain equilibrium 
as they cease to be relativistic. 

In 5 we discussed the quantum mechanics of a closely related problem including the relevent 
Fermi-Dirac and Einstein-Bose distributions. The Hamiltonian separates in the form H = H(5t) + 
H(r) + r~ 2 H(X.). The use of momenta allowed us to get the solutions without the introduction of 
r-time but the pulsations of the equilibrium have to be sought out in the correlations whereas they 
stand out more clearly in the classical case discussed above. 
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Figure 1: Variation of r 2 (top), total potential energy (middle) and V2 (bottom) during 10 pulsations 
following the first 50 pulsations. 
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Figure 2: Variation of the cluster shape E^/E 1 ?) a t the beginning of the run. Note that the 
initial anisotropy relaxes in about 5 pulsations 



7 



_l I I L_ 



10 



15 



t/per 



20 



Figure 3: Variation near the beginning of the run of the mean square peculiar velocities in the z(solid) 
and x(dashed) directions scaled with r. Note that these relax to equal values in about 5 pulsations 
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Figure 4: The equilibrium distribution of peculiar velocities. The upper graph shows the distribution 
of the peculiar velocities for eight phases of the pulsation cycle. The broad distributions correspond to 
phases when the cluster is compressed and the narrowest portions to expanded phases. In the middle 
graph the velocities have been rescaled with r demonstrating that the eight distributions coincide. The 
bottom portion demonstrates the Maxwellian nature of the distribution by showing the logarithm of the 
distribution function is linear in r 2 v 2 . 



